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Abstract. We discuss the statistical mechanics of violent relaxation in stellar sys- 
tems following the pioneering work of Lynden-Bell (1967). The solutions of the 
gravitational Vlasov-Poisson system develop finer and finer filaments so that a sta- 
tistical description is appropriate to smooth out the small-scales and describe the 
"coarse-grained" dynamics. In a coarse-grained sense, the system is expected to 
reach an equilibrium state of a Fermi-Dirac type within a few dynamical times. 
We describe in detail the equilibrium phase diagram and the nature of phase tran- 
sitions which occur in self-gravitating systems. Then, we introduce a small-scale 
parametrization of the Vlasov equation and propose a set of relaxation equations for 
the coarse-grained dynamics. These relaxation equations, of a generalized Fokker- 
Planck type, are derived from a Maximum Entropy Production Principle (MEPP). 
We make a link with the quasilinear theory of the Vlasov-Poisson system and derive 
a truncated model appropriate to collisionless systems subject to tidal forces. With 
the aid of this kinetic theory, we qualitatively discuss the concept of "incomplete 
relaxation" and the limitations of Lynden-BelPs theory. 



1 Introduction 

It has long been realized that galaxies, and self-gravitating systems in general, 
follow a kind of organization despite the diversity of their initial conditions 
and their environement. This organization is illustrated by morphological 
classification schemes such as the Hubble sequence and by simple rules which 
govern the structure of individual self-gravitating systems. For example, el- 
liptical galaxies display a quasi-universal luminosity profile described by de 
Vaucouleurs' i? 1 / 4 law and most of globular clusters are well fitted by the 
Michie-King model ||. The question that naturally emerges is, what deter- 
mines the particular configuration to which a self-gravitating system settles. 
It is possible that their present configuration crucially depends on the condi- 
tions that prevail at their birth and on the details of their evolution. However, 
in view of their apparent regularity, it is tempting to investigate whether 
their organization can be favoured by some fundamental physical principles 
like those of thermodynamics and statistical physics. We ask therefore if the 
actual states of self-gravitating systems are not simply more probable than 
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any other possible configuration, i.e. if they cannot be considered as max- 
imum entropy states. This thermodynamical approach may be particularly 
relevant for globular clusters and elliptical galaxies which are described by 
a distribution function that is almost isothermal || . In the case of globular 
clusters, the relaxation proceeds via two-body encounters and this collisional 
evolution is well-described by kinetic equations of a Fokker-Planck-Landau 
type for which a H-theorem is available. By contrast, for elliptical galaxies, 
two-body encounters are completely negligible (the corresponding relaxation 
time t co u exceeds the age of the universe by many orders of magnitude) 
and the galaxy dynamics is described by the Vlasov equation, i.e. collision- 
less Boltzmann equation Since the Vlasov equation rigorously conserves 
entropy, a relaxation towards an isothermal distribution looks at first sight 
relatively surprising. Yet, the inner regions of elliptical galaxies appear to be 
isothermal and this fact stemed as a mystery for a long time. 

In a seminal paper, Lynden-Bell argued that the violently changing 
gravitational field of a newly formed galaxy leads to a redistribution of ener- 
gies between stars and provides a mechanism analogous to a relaxation in a 
gas. The importance of this form of relaxation had previously been stressed 
by a number of authors including Henon and King but Lynden-Bell showed 
for the first time the relevance of a statistical description. He argued that the 
Vlasov-Poisson system develops an intricate mixing process in phase space 
associated with the heavily damped oscillations of a protogalaxy initially 
far from mechanical equilibrium and collapsing under its own gravity. As a 
result, the solutions of the Vlasov equation are not smooth but involve in- 
termingled filaments at smaller and smaller scales. In this sense, there is no 
convergence towards equilibrium but rather the formation of a fractal-like 
structure in phase space. However, if we introduce a macroscopic level of 
description and make a local average of the distribution function over the 
filaments, the resulting "coarse-grained" distribution function is smooth and 
is expected to reach a maximum entropy state (i.e. most mixed state) on a 
very short time scale of the order of the dynamical time to- This process is 
called violent relaxation and is acknowledged to account for the regularity of 
elliptical galaxies or other collisionless self-gravitating systems. Lynden-Bell 
predicted that the equilibrium state should be described by a Fermi-Dirac 
distribution function or a superposition of Fermi-Dirac distributions. Here, 
degeneracy is due not to quantum mechanics but to the Liouville theorem 
that prevents the smooth distribution function from exceeding the maximum 
of its initial value. In the non degenerate limit, the Fermi-Dirac distribution 
functions reduce to Maxwellians. The prediction of an isothermal distribution 



1 The H-theorem, proved by Boltzmann for an ideal gas, states that entropy is a 
monotonically increasing function of time. Boltzmann's kinetic theory of gases 
therefore provides a direct justification of the second principle of thermodynam- 
ics. The generalization of this theorem for self-gravitating systems rests on sim- 
plifying assumptions which are difficult to rigorously justify |3q] . 
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for collisionless stellar systems was considered as a triumph in the 1960's and 
laid the foundation of a new type of statistical mechanics. Of course, the va- 
lidity of the theory is conditioned by a hypothesis of ergodicity which may not 
be completely fulfilled. This is the complicated problem of "incomplete re- 
laxation" which limitates the power of prediction of Lynden-Bell's approach. 
However, as we shall see, these difficulties should not throw doubt on the 
importance of this statistical description. A similar relaxation process is at 
work in two-dimensional turbulence (described by the 2D Eulcr equation) 
and can explain the organization and maintenance of coherent vortices, such 
as the Great Red Spot of Jupiter, which are common features of large-scale 
geophysical or astrophysical flows |5^ , ^6| , ^3U5^ , ^6| ,|6|,[7| . The mathematical rel- 
evance of this statistical description has been given by Robert [ f50| introducing 
the concept of Young measures. The formal analogy between two-dimensional 
vortices and stellar systems has been discussed by Chavanis [pr|Jl^ , [l7| , p7| . 

This paper is organized as follows. In section 0, we introduce the gravi- 
tational Vlasov-Poisson system and list its main properties. In section ^, we 
present the statistical approach of Lynden-Bell 0| to the problem of violent 
relaxation. In section |], we show that the Fermi-Dirac equilibrium distri- 
bution predicted by Lynden-Bell is not entirely satisfactory since it has an 
infinite mass. We must therefore invoke incomplete relaxation and introduce 
truncated models. In the artificial situation in which the system is enclosed 
within a spherical box, we can calculate the Fermi-Dirac spheres explicitly 
and prove the existence of a global entropy maximum for each value of en- 
ergy. For low energies, this equilibrium state has a degenerate core surrounded 
by a dilute atmosphere, as calculated by Chavanis & Sommeria (25). More 
generally, we determine the complete equilibrium phase diagram and discuss 
the nature of phase transitions in self-gravitating systems. In section |^, we 
describe the coarse-grained relaxation of collisionless stellar systems towards 
statistical equilibrium in terms of a generalized Fokker-Planck equation. This 
relaxation equation is derived from a phenomenological Maximum Entropy 
Production Principle (MEPP) and involves a diffusion in velocity space com- 
pensated by a nonlinear friction. In section ^, we present a quasilinear theory 
of the Vlasov-Poisson system and show that it leads to a kinetic equation 
of a Landau type. When the system is close to equilibrium (so that a ther- 
mal bath approximation can be implemented) this equation reduces to the 
Fokker-Planck equation of the thermodynamical approach and the diffusion 
coefficient can be explicitly evaluated. This provides a new, self-consistent, 
equation for the "coarse-grained" dynamics of stellar systems where small 
scales have been smoothed-out in an optimal way. In section ^, we use this 
kinetic model to derive the distribution function of a tidally truncated col- 
lisionless stellar system. This truncated model preserves the main features 
of Lynden-Bell's distribution (including degeneracy) but has a finite mass, 
avoiding the artifice of a spherical container. Other truncated models at- 
tempting to take into account incomplete relaxation are discussed. 
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2 The gravitational Vlasov-Poisson system 

For most stellar systems, the encounters between stars are negligible @] and 
the galaxy dynamics is described by the self-consistent Vlasov-Poisson system 

df df df 

W + v ^ + F av-°' (1) 

A$ = 4ttG J /d 3 v. (2) 

Here, /(r, v, t) denotes the distribution function (defined such that fd 3 rd 3 v 
gives the total mass of stars with position r and velocity v at time t), 
F(r, t) = —V<P is the gravitational force (by unit of mass) experienced by 
a star and <£(r, t) is the gravitational potential related to the star density 
p(r,i) = J fd 3 v by the Newton-Poisson equation (||). The Vlasov equation 
(0) simply states that, in the absence of encounters, the distribution function 
/ is conserved by the flow in phase space. This can be written Df /Dt = 
where D/Dt = d/dt + XJqVq is the material derivative and U$ = (v, F) is a 
generalized velocity field in the 6-dimensional phase space (r, v) [by defini- 
tion, V6 = (d/dr,d/dv) is the generalized nabla operator]. Since the flow is 
incompressible, i.e. WqXJq — 0, the hypervolume of a "fluid" particle is con- 
served. Since, in addition, a fluid particle conserves the distribution function, 
this implies that the total mass (or hypervolume) of all phase elements with 
phase density between / and / + 6f is conserved. This is equivalent to the 
conservation of the Casimir integrals 



J h(f)d 3 vd 3 v, (3) 



for any continuous function h(f) (they include in particular the total mass 
M = J fd 3 rd 3 v). It is also straightforward to check that the Vlasov-Poisson 
system conserves the total energy (kinetic + potential) 

E = [ -fv 2 d 3 rd 3 v + - f f$d 3 rd 3 v = K + W, (4) 



T 2 
the angular momentum 



J / (r A v)d 3 rd 3 v, (5) 



and the impulse 



P = / fvd 3 rd 3 v. (6) 



In the following, we shall work in the barycentric frame of reference and 
assume that the system is non rotating so that the conservation of L and P 
can be ignored (see Refs. [Hp?]] for a generalization). 
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3 Lynden-Bell's approach of violent relaxation 

When the initial condition is far from equilibrium, the Vlasov-Poisson system 
develops a complicated mixing process in phase space and generates intermin- 
gled filaments due to stirring effects. Therefore, a deterministic description of 
the flow in phase space requires a rapidly increasing amount of information 
as time goes on. For that reason, it is appropriate to undertake a probabilistic 
description in order to smooth out the small-scales (fine-grained) and con- 
centrate on the locally averaged (coarse-grained) quantities. This statistical 
analysis has been considered a long time ago by Lynden-Bell Q. The sta- 
tistical mechanics of continuous systems is not as firmly established as in the 
usual case of N-body systems but Robert |5(J has recently developed a math- 
ematical justification of this procedure in terms of Young measures. We shall 
describe below the argumentation of Lynden-Bell, which is more intuitive in 
a first approach. 

Starting from some arbitrary initial condition, the distribution function is 
stirred in phase space but conserves its values 77 (levels of phase density) and 
the corresponding hypervolumes 7(77)^77 as a property of the Vlasov equation. 
Let us introduce the probability density p(r, v, 77) of finding the level of phase 
density 77 in a small neighborhood of the position r, v in phase space. This 
probability density can be viewed as the local area proportion occupied by 
the phase level 77 and it must satisfy at each point the normalization condition 

J p(r,v,?7)d?7 = 1. (7) 

The locally averaged distribution function is then expressed in terms of the 
probability density as 

/(r,v) = J p(r,v, 77)77^77, (8) 
and the associated (macroscopic) gravitational potential satisfies 

A$ = 4ttG [ /d 3 v. (9) 



Since the gravitational potential is expressed by space integrals of the density, 
it smoothes out the fluctuations of the distribution function, supposed at very 
fine scale, so <P has negligible fluctuations, ft is then possible to express the 
conserved quantities of the Vlasov equation as integrals of the macroscopic 
fields. These conserved quantities are the global probability distributions of 
phase density 7(77) (i.e., the total hypervolume occupied by each level 77 of 
phase density) 

7(77) = / p(r,v,77)d 3 rd 3 v, (10) 
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and the total energy 

E = J ^Jv 2 d 3 rd 3 v + i J J $d 3 rd 3 v. (11) 

As discussed above the gravitational potential can be considered as smooth, 
so we can express the energy in terms of the coarse-grained functions / and 
<P neglecting the internal energy of the fluctuations. 

To determine the equilibrium distribution of the system, we need to in- 
troduce an entropy functional like in ordinary statistical mechanics. As is 
customary, Lynden-Bell defines the mixing entropy as the logarithm of the 
number of microscopic configurations associated with the same macroscopic 
state (characterized by the probability density p(r,v,rf)). After introducing 
a counting "a la Boltzmann" , he arrives at the expression (4j| : 

S=-j p(r,v,ri)\np(r,v,r))drid 3 rd 3 v, (12) 

where the integral extends over phase space and over all the levels of phase 
elements. A mathematical justification of this entropy has been given by 
Robert . The most likely distribution to be reached at equilibrium is then 
obtained by maximizing the mixing entropy ( |l2| ) subject to the constraints 
([lO])([ll]) and the normalization condition (Q). This variational problem is 
treated by introducing Lagrange multipliers so that the first variations satisfy 



SS - f3SE - J a(ri)5-f(ri)dri - J C(r, v)5 \^ J pdr/j d 3 rd 3 v = 0, (13) 

where ft is the inverse temperature and a(r/) the "chemical potential" of 
species r\. The resulting optimal probability density is a Gibbs state which 
has the form 

e -a(r l )-f3r l (^+i>) 

v ' = - — , ; R jzr^ ■ (14) 

r e -a(»?)-/3>j(-V+ <p )d77 

The previous analysis gives a well defined procedure to compute the sta- 
tistical equilibrium states. The gravitational field is obtained by solving the 
Poisson equation (|^) with the distribution function (|J) determined by the 
Gibbs state (|T4|). The solution depends on the Lagrange multipliers (5 and 
a(r]) which must be related to the conserved quantities E and 7(77) by equa- 
tions (pi]) (10). This precedure determines critical points of entropy. Whether 



these critical points are maxima or not is decided by the sign of the second 
order variations of entropy. 

Following Lynden-Bell 0] , we consider a particular situation that presents 
interesting features and for which the previous problem can be studied in de- 
tail. Keeping only two levels / = 770 and / = is convenient to simplify the 
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discussion and is probably representative of more general cases. Within this 
"patch" approximation, the mixing entropy reduces to 

S =- //Il n I + (l_IWl_IU d 3 rd 3 V) (15) 

J Ivo Vo V W \ W J J 

and the equilibrium coarse-grained distribution / = p(r, v, T]q)t]o takes ex- 
plicitly the form 

7= 1 + A°e**" (16) 

where e = ^ + <P is the energy of a star (per unit of mass) and A = 
e a(vo)~a{o) > o an equivalent Lagrange multiplier. Equation Ji~6| ) is, apart 
from a reinterpretation of the constants, the distribution function for the 
self-gravitating Fermi-Dirac gas. Here, the exclusion principle / < T)q is due 
to the incompressibility constraint @) not to quantum mechanics. Because 
of the averaging procedure, the coarse-grained distribution function can only 
decrease by internal mixing, as vaccum is incorporated into the patch, and 
this results in an "effective" exclusion principle. Lynden-Bell's distribution 
therefore corresponds to a 4-th type of statistics since the "fluid" particles are 
distinghuishable but subject to an exclusion principle. However, formally, this 
distribution coincides with the Fermi-Dirac distribution. In the fully degen- 
erate case, this equilibrium has been studied extensively in connection with 
white dwarf stars Q . The structure of equilibrium may crucially depend on 
the degree of degeneracy (see section ||) but Lynden-Bell Q gives arguments 
according to which stellar systems would be non degenerate 0. In that limit, 
/ <C ijq, and equation ([l6|) reduces to the Maxwell-Boltzmann statistics 

J=Ae~ f)voe . (17) 

This was in fact the initial goal of Lynden-Bell in 1967: his theory of "vi- 
olent relaxation" was able to justify a Maxwell-Boltzmann equilibrium dis- 
tribution, without recourse to collisions, on a short time scale ~ to <C t C oii 
consistent with the age of ellipticals. In addition, the individual mass of the 
stars never appears in his theory based on the Vlasov equation. Therefore, 



the equilibrium state (17) does not lead to a segregation by mass in contrast 
with a collisional relaxation. This is in agreement with the observed light 
distributions in elliptical galaxies that would show greater colour differences 
if a marked segregation by mass was established. 



4 Computation of Fermi-Dirac spheres 

The first problem to be tackled is obviously the computation of self-gravitating 
Fermi-Dirac spheres. Let us first consider the non degenerate limit / <C rjo 

2 In fact, his arguments do not apply to galactic nuclei where his type of degeneracy 
may be important. 
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corresponding to a dilute system. In that case, the mixing entropy ( |l5| ) re- 
duces to the ordinary Boltzmann entropy 

/ln/d 3 rd 3 v, (18) 

whose maximization at fixed mass and energy leads to the Maxwcll-Boltzmann 
statistics ( |l7| ) . Substituting this optimal distribution in the Poisson equation 
(^|), we find that the gravitational potential is determined at equilibrium by 
the differential equation 

A$ = ^GA'e~^ , (19) 

where the Lagrange multipliers A' and f3 have to be related to the total mass 
and total energy of the system. The Boltzmann-Poisson equation ( |l9[) has 
been studied extensively in the context of isothermal gaseous spheres || and 
in the case of collisional stellar systems such as globular clusters (4f| . We can 
check by direct substitution that the distribution 

#.(r) = - H2nGpAr 2 ), p s {r) = (20) 

is an exact solution of equation ( |l9| ) known as the singular isothermal sphere 
[||. Since p ~ r -2 at large distances, the total mass of the system M = 
J^°° p4nr 2 dr is infinite! More generally, we can show that any solution of the 
Boltzmann-Poisson equation (|l9| ) behaves like the singular sphere as r — > +oo 
and has therefore an infinite mass. This means that no equilibrium can exist 
in an unbounded domain: the density can spread indefinitely while conserving 
energy and increasing entropy Q . 

In practice, the infinite mass problem does not arise if we realize that the 
relaxation of the system is necessarily incomplete [fl3|| . There are two major 
reasons for incomplete relaxation: (i) The mean field relaxation process is 
dependant on the strength of the variations in potential. As these die out, 
the relaxation ceases and it is likely that the system may find a stable steady 
state before the relaxation process is completed. Therefore, the relaxation is 
effective only in a finite region of space (roughly the main body of the galaxy) 
and during a finite period of time (while the galaxy is dynamically unsteady) . 
Orbits which lie partly outside the relaxing region and have periods longer 
than the time for which the galaxy is unsteady will not acquire their full 
quotas of stars, (ii) in practice, the galaxy will not be isolated but will be 
subject to the tides of other systems. Therefore, high energy stars will es- 
cape the system being ultimately captured by the gravity of a nearby object. 
These two independant effects have similar consequences and will produce 
a modification of the distribution function at high energies. Some truncated 
models can be obtained by developing a kinetic theory of encounterless relax- 
ation (sections pHfl). This kinetic approach will provide a precise framework 
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to understand what limits relaxation and why complete equilibrium is not 
reached in general. 

However, for the present, we shall avoid the infinite mass problem by 
confining artificially the system within a box of radius R. The calculation of 
finite isothermal spheres has been carried out by Antonov EJ|, Lynden-Bell 
& Wood pi, Katz ||, Padmanabhan de Vega & Sanchez @ and 

Chavanis |L8| (see Chavanis [^l) for an extension in general relativity) . These 
studies were performed in the framework of gaseous stars or collisional stellar 
systems (e.g., globular clusters) but they extend in principle to violently 
relaxed collisionlcss stellar systems described by the distribution function 
([l7|). Therefore, we shall give a brief summary of these classical results before 
considering the case of the Fermi-Dirac distribution (|l^) . The phase diagram 
is represented in Fig. [l] where the inverse temperature is plotted as a function 
of minus the energy. It is possible to prove the following results: (i) there 
is no global maximum for the Boltzmann entropy, (ii) there are not even 
critical points for the Boltzmann entropy if A = —ER/GM 2 > 0.335. (iii) 
local entropy maxima (LEM) exist if A = -ER/GM 2 < 0.335; they have a 
density contrast 1Z — p(0)/p(R) < 709 (upper branch), (iv) critical points 
of entropy with density contrast 1Z > 709 are unstable saddle points SP 
(spiraling curve). Conclusions (i) and (ii) have been called "gravothermal 
catastrophe" or "Antonov instability". When A = -ER/GM 2 > 0.335, there 
is no hydrostatic equilibrium and the system is expected to collapse and 
overheat. This is a natural evolution (in a thermodynamical sense) because 
a self-gravitating system can always increase entropy by taking a "core-halo" 
structure and by making its core denser and denser (and hotter and hotter) 
[|J. As discussed by Lynden-Bell & Wood |l5|, this instability is probably 
related to the negative specific heats of self-gravitating systems: by losing 
heat, the core grows hotter and evolves away from equilibrium. 

The "gravothermal catastrophe" picture has been confirmed by sophisti- 
cated numerical simulations of globular clusters that introduce a precise de- 
scription of heat transfers between the "core" and the "halo" using moment 
equations 0] , orbit averaged Fokker-Planck equation j29| or fluid equations 
for a thermally conducting gas [Q. In these studies, the collapse proceeds 
self-similarly (with power law behaviours) and the central density becomes 
infinite in a finite time. This singularity has been known as "core collapse" 
and many globular clusters have probably experienced core collapse. This is 
certainly the most exciting theoretical aspect of the collisional evolution of 
stellar systems. In practice, the formation of hard binaries can release suffi- 
cient energy to stop the collapse and even drive a reexpansion of the system 
[ p2| . Then, a series of gravothermal oscillations should follow but it takes 
times much larger than the age of globular clusters ||. It has to be noted 
that, although the central density tends to infinity, the mass contained in the 
core tends to zero so, in this sense, the gravothermal catastrophe is a rather 
unspectacular process S. However, the gravothermal catastrophe may also 
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Fig. 1. Equilibrium phase diagram for classical isothermal spheres. 



be at work in dense clusters of compact stars (neutron stars or stellar mass 
black holes) embedded in evolved galactic nuclei and it can presumably initi- 
ate the formation of supermassive black holes via the collapse of such clusters. 
Indeed, when the central redshift (related to the density contrast) exceeds a 
critical value z c ~ 0.5, a dynamical instability of general-relativistic origin 
sets in and the star cluster undergoes a catastrophic collapse to a black hole 
on a dynamical time |Q ■ It has been suggested that this process leads nat- 
urally to the birth of supermassive black holes of the right size to explain 



quasars and AGNs |56 



Since collisionless stellar systems undergoing violent relaxation are de- 
scribed by an isothermal distribution function of the Fermi-Dirac type (|l6|), 
it is particularly relevant to ask whether galaxies in their birth stages can 
undergo a form of gravothermal catastrophe and if so what they do about it. 
The collapse of galaxies was described qualitatively by Lynden-Bell & Wood 
p5| : "...The centers will then begin to separate into a core - a sort of separate 
body which might even be called a nucleus. This will cease to shrink when 
it becomes degenerate in Lynden-Bell's sense. The system will then have a 
core-halo structure which is an equilibrium of an isothermal Fermi-Dirac gas 
sphere. These will show a variety of different nuclear concentrations depend- 
ing on the degeneracy parameter. It is evident that theory developed along 
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these lines has the chance of making sense of the variety of different galactic 
nuclei." 
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The computation of isothermal Fermi-Dirac spheres has been performed 
only recently by Chavanis & Sommeria |p5[ . It can be proved that a global 
entropy maximum exists for all accessible values of energy |E5 51 . Therefore, 



degeneracy has a stabilizing role and is able to stop the "gravothermal catas- 
trophe" . The equilibrium diagram is represented i n Figs. and now depends 
on the degeneracy parameter /i = 770/ (/) = ??oV '5127r 4 G 3 M R 3 , where 770 is 
the maximum allowable value of the distribution function and (/) its typical 
average value in the box of radius R. We see that the inclusion of degeneracy 
has the effect of unwinding the spiral (dashed line). When /1 is small (Fig. 
^|) , there is only one critical point of entropy for each value of energy and it 
is a global entropy maximum (GEM). For small values of A (high energies) 
the solutions almost coincide with the classical, non degenerate, isothermal 
spheres. When A is increased (low energies) the solutions take a "core-halo" 
structure with a partially degenerate core surrounded by a dilute Maxwcllian 
atmosphere. It is now possible to overcome the critical energy A c — 0.335 
and the critical density contrast 1Z ~ 709 found by Antonov for a classical 
isothermal gas. In that region, the specific heat C — dE/dT is negative. As 
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energy decreases further, more and more mass is concentrated into the nucleus 
(which becomes more and more degenerate) until a minimum accessible en- 
ergy, corresponding to A max (fj,), at which the nucleus contains all the mass. In 
that case, the atmosphere has been "swallowed" and the system has the same 
structure as a cold white dwarf star || . This is a relatively singular limit since 
the density drops to zero at a finite radius whereas for partially degenerate 
systems, the density decays like r~ 2 at large distances. When the degeneracy 
parameter fj, is large (Fig. [|), there are now several critical points of entropy 
for each single value of energy in the range /t*(/i) < A < A c . The solutions on 
the upper branch (points A) are non degenerate and have a smooth density 
profile; they form the "gaseous" phase. The solutions on the lower branch 
(points C) have a "core-halo" structure with a degenerate nucleus and a di- 
lute atmosphere; they form the "condensed" phase. According to the theorem 
of Katz |39), they are both entropy maxima while the intermediate solutions, 
points B, are unstable saddle points (SP). These points are similar to points 
A, except that they contain a small embryonic nucleus (with small mass and 
energy) which plays the role of a "germ" in the langage of phase transitions. 
For small values of A, the non degenerate solutions (upper branch) are global 
entropy maxima (GEM). They become local entropy maxima (metastable 
states) at a certain A t > A* depending on the degeneracy parameter. At 
that A t , the degenerate solutions (lower branch) that were only local entropy 
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maxima (LEM) suddenly become global entropy maxima. We expect there- 
fore a phase transition to occur between the "gaseous" phase (upper branch) 
and the "condensed" phase (lower branch) at A t . For A t ([i) < A < A c , the 
non degenerate solutions (points A) are metastable but we may suspect [^5| 
that they are long-lived so that they are physical. It is plausible that these 
metastable states will be selected by the dynamics even if states with higher 
entropy exist. This depends on a complicated notion of "basin of attraction" 
as studied by Chavanis, Rosier and Sire |22| with the aid of a simple model 
of gravitational dynamics. Therefore, the true phase transition will occur at 
A c = 0.335 at which the gaseous phase completely disappears: at that point, 
the gravothermal catastrophe sets in but, for collisionless stellar systems (or 
for fermions), the core ultimately ceases to shrink when it becomes degen- 
erate. In that case, the system falls on to a global entropy maximum which 
is the true equilibrium state for these systems (point D). This global en- 
tropy maximum has a "core-halo" structure with a degenerate core and a 
dilute atmosphere. A simple analytical model of these phase transitions has 
been proposed in Ref. 1 20 and provides a fairly good agreement with the full 
numerical solution. A particularity of self-gravitating systems, which are in 
essence non-extensive, is that the statistical ensembles (microcanonical and 
canonical) are not interchangeable. Therefore, the description of the equilib- 
rium diagram is different whether the system evolves at fixed energy of fixed 
temperature. A discussion of this interesting phenomenon can be found in 
the review of Padmanabhan [till and in Chavanis [E0| . 



For astrophysical purposes, it is still a matter of debate to decide whether 
collisionless stellar systems like elliptical galaxies are degenerate (in the sense 
of Lyndon-Bell) or not. Since degeneracy can stabilize the system without 
changing its overall structure at large distances, we have suggested that de- 
generacy could play a role in galactic nuclei The recent simulations of 
Leeuwin and Athanassoula (42) and the theoretical model of Stiavelli (57J are 
consistent with this idea especially if the nucleus of elliptical galaxies contains 
a primordial massive black hole. Indeed, the effect of degeneracy on the dis- 
tribution of stars surrounding the black hole can explain the cusps observed 
at the center of galaxies |37|. This form of degeneracy is also relevant for 
massive neutrinos in Dark Matter models where it competes with quantum 
degeneracy [^0|. In fact, the Fermi-Dirac distribution of massive neutrinos 
in Dark Matter models (which form a collisionless self-gravitating system) 
might be justified more by the process of "violent relaxation" |43}] than by 
quantum mechanics [|34"p|] . As shown by Chavanis & Sommeria ]25[ , violent 
relaxation can lead to the formation of a dense degenerate nucleus with a 
small radius and a large mass. This massive degenerate nucleus could be an 
alternative to black holes at the center of galaxies (|flE5| . 
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5 The Maximum Entropy Production Principle 

Basically, a collisionless stellar system is described in a self-consistent mean 
field approximation by the Vlasov-Poisson system In principle, these 

coupled equations determine completely the evolution of the distribution 
function /(r, v,t). However, as discussed in section ^, we are not interested 
in practice by the finely striated structure of the flow in phase space but only 
by its macroscopic, i.e. smoothed-out, structure. Indeed, the observations 
and the numerical simulations are always realized with a finite resolution. 
Moreover, the "coarse-grained" distribution function / is likely to converge 
towards an equilibrium state contrary to the exact distribution / which de- 
velops smaller and smaller scales. 

If we decompose the distribution function and the gravitational potential 
in a mean and fluctuating part (J = f + f, & = <P + and take the local 
average of the Vlasov equation (|lj), we readily obtain an equation of the form 

dj df -dj dJ f 

for the "coarse-grained" distribution function with a diffusion current J/ = 

/F related to the correlations of the "fine-grained" fluctuations. 

The problem in hand consists in determining a closed form for the dif- 
fusion current J/. Its detailed expression fF depends on the subdynamics 
and is therefore extremely difficult to capture from first principles. Indeed, 
the "violent relaxation" is a very nonlinear and very chaotic process rul- 
ing out any attempt to implement perturbative methods if we are far from 
equilibrium. For that reason, there is a strong incentive to explore varia- 
tional methods which are considerably simpler and give a more intuitive un- 
derstanding of the problem. We propose to describe the relaxation of the 
coarse-grained distribution function / towards the Gibbs state (|l4| ) with a 
Maximum Entropy Production Principle [p7[ . This thermodynamical prin- 
ciple assumes that: "during its evolution, the system tends to maximize its 
rate of entropy production S while satisfying all the constraints imposed by 
the dynamics" . There is no rigorous justification for this principle and it is 
important therefore to confront the MEPP with kinetic theories (when they 
are available) in order to determine its range of validity. In any case, the 
MEPP can be considered as a convenient tool to build relaxation equations 
which are mathematically well-behaved and which can serve as numerical 
algorithms to calculate maximum entropy states. It is also "perfectible" in 
the sense that any new information about the dynamics of the system can 
be incorporated. This principle is reminiscent of Jayncs Q ideas and is a 
clear extension of the well-known principle of statistical mechanics according 
to which: "at equilibrium, the system is in a maximum entropy state con- 
sistent with all the constraints" . This is also the most natural extension of 
Lynden-Beh"s theory out of equilibrium. 
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Like for equilibrium states, the evolution of the flow in phase space is 
described by a set of local probabilities p(r,v,rj,t) and the locally averaged 
gravitational potential ^(r, t) is obtained by the integration of equation (||) 
where f(r,v,t) is replaced by /(r,v,i) = J p(r,v,T],t)r]dr]. The phase ele- 
ments are thus transported in phase space by the corresponding averaged 
velocity field Ve = (v, — V<P) and we suppose that, in addition, they undergo 
a diffusion process. This diffusion occurs only in velocity space, due to the 
fluctuations of the gravitational field There is no diffusion in position space 
since the velocity v is a pure coordinate and does not fluctuate. As a result, 
the probability densities will satisfy a convection-diffusion equation of the 
form 

dp ^ dp ^ — dp dJ 

dt dr 9v dv' 
where J(r, v, 77, t) is the diffusion current of the phase element r\. This equa- 
tion conserves the total hypervolumc (in phase space) occupied by each phase 
clement. Multiplying equation (22) by r\ and integrating over the levels of 
phase density returns equation (21) with Jj = J J(r,v,r],t)f]di]. 

To apply the MEPP, we need first to compute the rate of change of en- 
tropy during the convection-diffusion process. Taking the time derivative of 
equation jl^), expressing dtp by equation ( [22] ) and noting that phip is con- 
served by the advective term, we get 

S = - J J^d 3 Td 3 vd V . (23) 

We now determine J such that, for a given field p at each time t, J maximizes 
the entropy production S with appropriate dynamical constraints, which are: 
-the conservation of the local normalization condition (Q) implying 

fj(r,v, 77,^77 = 0, (24) 

-the conservation of energy expressed from equations (11) and ( ^l|) as 



E = J J/V<2 3 rcZ 3 v = 0, (25) 

- a limitation on the eddy flux |J|, characterized by a bound C(r, v, t), which 
exists but is not specified 

^-*7<C(r,v,*). (26) 
2p 

This variational problem can be solved by introducing (at each time t) La- 
grange multipliers f, /?, 1/D for the three respective constraints. It can be 
shown by a convexity argument that reaching the bound ( |26| ) is always favor- 
able for increasing S, so that this constraint can be replaced by an equality. 
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The condition 



5S-135E- J jj S (^J ^-e^d 3 rd 3 v- f t 6 (J ^dr^d 3 rd 3 v = 0, (27) 



yields an optimal current of the form 



-D(r,v,t) 



dp 

dv 



(28) 



The Lagrange multiplier £ has been eliminated, using the condition ( p4[ ) of 
local normalization conservation. At equilibrium, the diffusion currents must 
vanish and we can check that this yields the Gibbs state ( |l4| ) with /3 the 
corresponding inverse temperature During the evolution, this quantity 
varies with time and is determined by the condition of energy conservation 
IBI). Introducing ( |2§| ) in the constraint (plf), we get 



fD^vd 3 rd 3 v 
J D{p - J 2 )v 2 d 3 rd 3 v' 



(29) 



where /" = J prf'drj are the local moments of the fine-grained distribution 
function. We have thus obtained a complete set of relaxation equations which 
exactly conserve the distribution of phase levels and energy. The rate of 
entropy production can be put in the form p7f : 



/2 

S = I — d 3 vd 3 vdr 1 , 
Dp 



(30) 



so the diffusion coefficient D must be positive for entropy increase. Except for 
its sign, the diffusion coefficient is not determined by the thermodynamical 
approach as it is related to the unknown bound C(r, v, t). 

These relaxation equations (|2^) (|2S|) (P9|) can be simplified in the case of 
a single non zero density level 770 and provide a new non linear equation for 
the evolution of the coarse-grained distribution function / = prjo : 



df df df d 
— + v— + F— = — 
dt dr dv dv 



£>(^ + /3/(%-/)v 



(31) 



In the non degenerate limit (/ <C 770), equation 
Fokker-Planck equation 



51) takes the form of a 



df df df d 
— + v— + F— = — 
dt dr dv dv 



(32) 



This Fokker-Planck equation (sometimes called the Kramers-Chandrasekhar 
equation) is well-known in the case of collisional stellar systems (without the 
bar on /!) and is usually derived from a Markov hypothesis and a stochastic 
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Langevin equation [ fL0| . It can also be obtained from the following argument: 
because of close encounters, the stars undergo brownian motion and diffuse 
in velocity space (responsible for the first term in the right hand side) . How- 
ever, under the influence of this diffusion, the kinetic energy per star will 
diverge as < v 2 /2 >~ 3Dt and one is forced to introduce an "ad hoc" dy- 
namical friction (second term in the right hand side) with a friction coeffi- 
cient £ = PDrjo (Einstein relation) in order to compensate this divergence 
and recover the Maxwellian distribution of velocities at equilibrium. From 
these two considerations (which express the fluctuation-dissipation theorem) 
results the ordinary Fokker- Planck equation (|32|). The Fokker-Planck equa- 
tion is justified here by an argument of a very wide scope that does not 
directly refer to the subdynamics: it appears to maximize the rate of entropy 
production at each time with appropriate constraints. Therefore, it applies 
equally well to collisional or (coarse-grained) collisionless stellar systems. In 
this description, the diffusion term and the "dynamical friction" directly re- 
sult from the variational principle and are associated with the variations of 
S and E respectively. Furthermore, the inverse temperature (3 appears as a 
Lagrange multiplier associated with the conservation of energy and the Ein- 
stein relation is automatically satisfied. In addition, our procedure can take 
into account the degeneracy of collisionless stellar systems, keeping equation 
( |3l| ) instead of the non degenerate limit (|32|). In the degenerate case, the fric- 
tion is non linear in / so that equation ( pl[ ) is not a Fokker-Planck equation 
in a strict sense. This nonlinearity is necessary to recover the Fermi-Dirac 
distribution function at equilibrium (while a linear friction drives the system 
towards the Maxwell-Boltzmann distribution). This generalization is impor- 
tant because degeneracy is specific to collisionless systems and may be crucial 
for the existence of an equilibrium state (see section ^) . 

In equations ([H]) (|3^), (3 is not constant but evolves with time so as to 
conserve the total energy. In the non degenerate limit, and assuming that D 
is constant, we get after a part integration 

f #£vd 3 rd 3 v 3M , . 

(3 t = - % 2M „ = — — 33 
J T)ofv 2 (PrcPv 2r] K{t) 

This equation relates the formal Lagrange multiplier (3{t) to the inverse of 
the average kinetic energy. This is of course satisfying on a physical point 
of view. An alternative Fokker-Planck equation involving a local temperature 
T(r, t) instead of /3(i) _1 has been proposed by Clemmow & Dougherty [^8| in 
the case of collisional systems. The energy is assumed to be locally conserved 
by the collisions, which is valid when the mean free path is much smaller 
than the size of the system. By contrast, this hypothesis does not seem to be 
justified for the violent relaxation of a collisionless system, which is rather a 
global process. 

Let us briefly review the main properties of our relaxation equations. First 
of all, they rigorously satisfy the conservation of energy and phase space hy- 
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pervolumes like the Vlasov-Poisson system (the conservation of impulse and 
angular momentum can also be satisfied by introducing appropriate Lagrange 
multipliers ^7j). Moreover, they guarantee the increase of entropy at each 
time (S > 0) with an optimal rate. Of course, this H-theorem is true for the 
coarse-grained entropy S = — J phi pd 3 rd 3 vdr) and not for the fine-grained 
entropy Sf. g = — J /ln/d 3 rd 3 v which is constant, as the integral of any 
function of /. The source of irreversibility is due to the coarse-graining proce- 
dure that smoothes out the small-scales and erases the microscopic details of 
the evolution. Accordingly, the relaxation equations (f22f) (^8|) ( |2£ ) a re likely to 



drive the system towards an equilibrium state (the Gibbs state (14)) contrary 
to the Vlasov equation that develops finer and finer scales. In mathematical 
terms, this means that the distribution function / converges to an equilib- 
rium state f in a weak sense. In fact, the situation is more complicated since 
the Gibbs state does not exist in an infinite domain (section ||). We shall 
see, however, in section [?] that the diffusion coefficient is proportional to the 
fluctuations of the distribution function (equation (|69|)). Now, in physical 
situations, these fluctuations vanish before the systems had time to relax 
completely. As a result, the relaxation stops and the system remains frozen 
in a subdomain of phase space. It is only in this subdomain (corresponding to 
the main body of the galaxy) that the Gibbs state is justified. This provides 
a physical mechanism for confining galaxies and justifying truncated mod- 
els. Alternatively, if the galaxy is not isolated but subject to the tides of a 
neighboring object, a tidally truncated model can be explicitly derived from 
these relaxation equations (section ^). It has a finite mass while preserving 
the essential features of Lynden-BelPs distribution function. 



The relaxation equations (122) (Eq)(29) can be simplified in two particular 



limits when: (i) the initial condition is approximated by a single level of phase 
density 770 surrounded by vaccum (ii) degeneracy, in Lyndon-Bell's sense, is 
neglected. These simplifications lead to the Fokker-Planck equation (|3^). If 
the galaxy has sufficiently large energy, this equation will drive the system 
towards a Maxwell-Boltzmann equilibrium state (|l7|). However, if the energy 
falls below a critical value, the Fokker-Planck equation does not reach any 
equilibrium state anymore and the system can achieve ever increasing values 
of entropy by developing core collapse (Antonov instability). For collision- 
less stellar systems, this "gravothermal catastrophe" should stop when the 
center of the galaxy becomes degenerate (see section ||). In that case, the 
Fokker-Planck equation ([32]) is not valid anymore and must be replaced by 
the degenerate relaxation equation ([n]). This equation converges towards the 
Fermi-Dirac equilibrium state (|l^), which exists for all values of energy. 

In summary, the MEPP is able to yield relevant relaxation equations 
for the coarse-grained dynamics of collisionless stellar systems experiencing 
violent relaxation. This relatively elegant and simple variational principle 
shows that the global structure of the relaxation equations is determined by 
purely thermodynamical arguments. All explicit reference to the subdynamics 
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is encapsulated in the diffusion coefficient which cannot be captured by the 
MEPP (it appears as a Lagrange multiplier related to an unknown bound on 
the diffusion flux) . ft must be therefore calculated with a kinetic model such 
as the quasilinear theory described in the next section. 



6 The quasilinear theory 

The quasilinear theory of the Vlasov-Poisson system was first considered 
by Kadomtsev & Pogutse |37|] fo r a homogeneous Coulombian plasma and 
extended by Severne & Luwel jp5[ for an inhomogeneous gravitational system. 
This theory was further discussed by Chavanis Jl3| in an attempt to make a 
link with the Maximum Entropy Production Principle. We shall give here a 
simple account of the quasilinear theory. Further details can be found in Refs. 
[ ^7|j55| , p^t . Our objective is to obtain an expression for the diffusion current 
J by working directly on the Vlasov-Poisson system, i.e. without assuming 
that the entropy increases as is done in the thermodynamical approach. 

Since the diffusion current J/ = /F is related to the fine-grained fluctu- 
ations of the distribution function, any systematic calculation starting from 
the Vlasov equation (|]J) must necessarily introduce an evolution equation 
for /. This equation is simply obtained by substracting equation ([H]) from 
equation (pi). This yields 



df df -df ~df ~df ^df 

— + v— + F — = -F— - F— + F— . (34) 

dt dr dv dv dv dv 

The essence of the quasilinear theory is to assume that the fluctuations are 
weak and neglect the nonlinear terms in equation ([}4|) altogether. In that 
case, equations pi"| ) and (34) reduce to the coupled system 



df - d — 

■s + (35) 

f < 36 > 

where L = + F^ is the advection operator in phase space. Physically, 
these equations describe the coupling between a subdynamics (here the small 
scale fluctuations /) and a macrodynamics (described by the coarse-grained 
distribution function /). 
Introducing the Greenian 

G{t 2 , h) = expJ - / dtL(t) \ , (37) 



we can immediately write down a formal solution of equation (|36|) , namely 
/(r,v,t)=G(f,0)/(r,v,0)-y dsG(t, t - s)F(r, t- s)^(r, v, t - s). (38) 
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Although very compact, this formal expression is in fact extremely compli- 
cated. Indeed, all the difficulty is encapsulated in the Greenian G(t, t — s) 
which supposes that we can solve the smoothed-out Lagrangian flow 

between t and t — s. In practice, this is impossible and we will have to make 
some approximations. 

The objective now is to substitute the formal result J38| ) back into equa- 
tion ( |35| ) and make a closure approximation in order to obtain a self-consistant 
equation for /(r, v, t). If the fluctuating force F were external to the system, 
we would simply obtain a diffusion equation 

%+Lj=^-(D^f-), (40) 
with a diffusion coefficient given by a Kubo formula 

D» v = [ dsF»(r,t)F»(r,t- s). (41) 
Jo 

However, in the case of the Vlasov-Poisson system, the gravitational force is 
in fact produced by the distribution of matter itself and this coupling will 
give rise to a friction term in addition to the pure diffusion. Indeed, we have 

F(r, t) = J F(r' -> r)/(r', v', t)dVdV, (42) 



where 



F ( r '- r )=G^4, (43) 



represents the force created by a (field) star in r' on a (test) star in r (New- 
ton's law). Therefore, considering equations ( |38| ) and (ff2|), we see that the 
fluctuations of the distribution function /(r, v, t) are given by an iterative 
process: f(t) depends on F(i — s) which itself depends on f(t — s) etc... We 
shall solve this problem perturbatively in an expansion in powers of the grav- 
itational constant G. This is the equivalent of the "weak coupling approxima- 
tion" in plasma physics. To order G 2 , we obtain after some rearrangements 



dt 



x \f v (v" r)/(r', v',t- s)f(r", v",t- s)^;{v, v,t-s) 

df 



+F v {v" - r')/(r, v, t - s)f(r", v", t - s)^j(r', v', t - s) 



(44) 
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In this expression, the Greenian G refers to the fluid particle r(t), v(£) and the 
Greenian G' to the fluid particle r'(t), v'(t). To close the system, it remains for 
one to evaluate the correlation function /(r, v, i)/(r', v', t). We shall assume 
that the mixing in phase space is sufficiently efficient so that the scale of the 
kinematic correlations is small with respect to the coarse-graining mesh size. 
In that case, 

/(r,v,i)/(r',v',i) = e 3 r e 3 v 6(r - r')S(v - v')p(r, v, t), (45) 

where e r and e v are the resolution scales in position and velocity respectively. 
Now, 

J 2 = (I -1? =T -f ■ (46) 

We shall assume for simplicity that the initial condition in phase space con- 
sists of a patch of uniform distribution function (/ = 770) surrounded by 
vaccum (/ = 0). This is the two-levels approximation already considered in 
sections || and |[ In that case f 2 — r/Q x / = r/ f and, therefore, 

/(r,v, *)/>', v',t) = 445(r - r')5(v - V)J{ m - J). (47) 

Substituting this expression in equation ( [i~4| ) and carrying out the integrations 
on r" and v", we obtain 

^ + Lj= eUl-^ J ds J d 3 r'd 3 v'F»(r' - v) t F»(r' r) t _ s 

X {7(Va-7)§;-7(m-J)§t} t _- (48) 

We have written j' t _ s = J(r'(t - s), v'(f - s),t-s), J t _ s = J(r(t - s), v(t - 
s),t-s), F"(r' -» t) t = F»(r'(t) -» r(t)) and i^(r' -» r) t _ s = F v (r'(t-s) -> 
r(t — s)) where r(t — s) and v(t — s) are the position and velocity at time 
t — s of the stellar fluid particle located in r = r(t), v = v(t) at time It is 
determined by the characteristics (^) of the smoothed-out Lagrangian flow. 

Equation ( [i"8| ) is a non Markovian integrodifferential equation: the value 
of / in r, v at time t depends on the value of the whole field f(r', v', t — s) 
at earlier times. If the decorrelation time t is short, we can make a Markov 
approximation and replace the bracket at time t — s by its value taken at time 
t. Noting furthermore that the integral is dominated by the contribution of 
field stars close to the star under consideration (i.e. when r' — > r), we shall 

make a local approximation and replace / (770 — / ) and by their values 
taken at r. In that case, the foregoing equation simplifies in 

H + Lf = ^ J dWv'F'V -> r) t FV - r) t _ s 

x{f'(Vo-f')§;-f(Vo-7)§^), (49) 
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where, now, / = /(r,v',t). The explicit reference to the past evolution of 
the system is only retained in the memory function 



J ds J dVF"(r' r) t F v {r' r) t 



This function can be calculated explicitly if we assume that, between t — s and 
t, the stars follow linear trajectories, so that v(t— s) = v and r(t— s) = r — vs 
[||. This leads to the generalized Landau equation 

f + i7 - 4? I "«~{r<m -r>§- n» -7)^}, (») 

where if '"' is the tensor 

Jf"" = 2vG 2 e 3 r e 3 v hi A- ( 6» v - ^f- ), (51) 
u \ u 1 ) 

and u = v' — v, ln/1 = ln(i?/e r ). This equation applies to inhomogeneous 
systems but, as a result of the local approximation, the effect of inhomogene- 
ity is only retained in the advective term. Equation (|5^) is very similar to the 
well-known Landau equation of collisional self-gravitating systems and elec- 
tric charges 38 . There are nevertheless two important differences: (i) The 
friction term is non linear in /, accounting for the degeneracy discovered by 
Lynden-Bell at equilibrium, (ii) The diffusion coefficient is proportional to 
the mass r/oe^e^, of a macrocell completely filled by the phase fluid instead of 
the mass m of an individual star in the ordinary Landau equation. In general 
77o£^e^ >mso that the relaxation by phase mixing is much more rapid than 
the collisional relaxation. From the above theory, we find that the collision- 
less relaxation operates on a time scale ~ tjj, the dynamical time, whereas 
the collisional relaxation operates on a time scale t co u ~ j^jf^D 3> tn where 
N ~ 10 12 is the number of stars in the galaxy. Therefore, the relaxation by 
phase mixing really corresponds to a "violent relaxation" f43| ] . 

It can be shown that the generalized Landau equation (|50|) conserves mass 
and energy 0. In fact, as a result of the local approximation, the energy is 
conserved locally, as if the system were homogeneous, and we have 

§0 y d3v = / J/vrf3v = °' Vr - (52) 

It is also easy to show that equation ( |50| ) satisfies a H-theorem (S > 0) for the 
Fermi-Dirac entropy (]l5|) . When a stationary state is reached S = and the 
Fermi-Dirac distribution ( fl6| ) is obtained, in agreement with Lynden-BelPs 
statistical theory p3| ]. This provides therefore another way of justifying his 
results from a dynamical point of view which does not explicitly rely on a 
maximization of entropy (the 77-theorem is not assumed but derived from 
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the kinetic theory). It can be noted that these properties result from the 
symmetry of the Landau collision term and not from Lagrange multipliers 
like in the thermodynamical approach. This is more satisfactory on a physical 
point of view. 

It is important to stress, however, that this quasilinear theory cannot 
describe the early, very non linear and very chaotic, stages of the "violent 
relaxation" . Indeed, the detailed study of Severne & Luwel ]55| reveals that 
the various approximations introduced in the theory make equation (50) ap- 
plicable only for t 3> tjj, where tjj is the dynamical time. Since the relaxation 
time is precisely of order to, the quasilinear theory is only marginally appli- 
cable and can describe, at most, the late quiescent phases of the relaxation, 
when the fluctuations have weaken. By contrast, it is plausible that equation 
( |3l| ) is more general and more appropriate to the context of "violent relax- 
ation" since it results from a thermodynamical approach which exploits at 
best the chaoticity of the system and the complete lack of information that 
we have to face at small scales. As a clear difference, it should be noted that 
the MEPP takes into account only the global conservation of energy (the La- 
grange multiplier j3{t) is a kind of inverse average temperature determined 
by an integration over the whole system) whereas the approximations intro- 
duced in the kinetic model lead to a local conservation of energy. This strong 
locality cannot account for the rather collective processes which are involved 
in the violent relaxation and may unveil a failure of the quasilinear theory. 
However, in section [7], we show that the two approaches are consistant if 
we are close to equilibrium and we use the quasilinear theory as a model to 
determine an explicit expression for the diffusion coefficient that appears in 
equation (pi 



7 Truncated models for collisionless stellar systems 

The equations provided by the MEPP and by the quasilinear theory have 
a very different mathematical structure. The relaxation equation (J3l]) is a 
partial differential equation whereas the generalized Landau equation ( |50| ) 
is an integrodifferential equation: the value of /(v) at time t + dt depends 
explicitly on the whole distribution of velocities /(v') at time t through an 
integration over v'. The usual way to transform an integrodifferential equa- 
tion into a partial differential equation is to make a guess for the function 
appearing in the integral and refine the initial guess by successive iterations. 
In practice, we simply make one sensible guess. Therefore, if we are close to 
equilibrium (and this is in fact dictated by the conditions of validity of the 
quasilinear theory), it seems natural to replace the distribution function / 
by the Fermi-Dirac distribution 

r = — 9° (53) 

l + \ e 0no — 
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In more physical terms, this amounts to a "thermal bath approximation" : the 
stars have not yet relaxed completely, but when we focus on the relaxation 
of a given stellar fluid particle (described by /) we can consider, in a first 
approximation, that the rest of the system (described by / ) is at equilibrium. 
With this thermal bath approximation, the generalized Landau equation ( j50| ) 
reduces to the nonlinear partial differential equation Jl3|: 



df df -df d 
— + v— + F— = — 
dt dr dv <9v 



(54) 



with a diffusion coefficient 



Equation ( |54| ) is precisely the equation derived from the MEPP. Together 
with the explicit expression of the diffusion coefficient ( pq ) it provides a self- 
consistent model for the "coarse-grained" dynamics of collisionless stellar 
systems experiencing violent relaxation. More general equations including 
anisotropic effects can also be obtained from this formalism jl3| . 

We shall now assume that the galaxy is not isolated but subject to the 
tides of other systems. In that case, high energy stars that have elongated 
orbits are removed by the gravity of these objects. We seek therefore a station- 
ary solution of equation (|54| ) satisfying the boundary condition /(e m ) = 0, 
where e m is the escape energy above which / = 0. This solution will provide 
a truncated distribution function with a finite mass |T^| . In fact, this problem 
was already tackled by Lynden-Bell in his seminal paper f4l| . He proposed 
to describe the evolution of the coarse-grained distribution function / by the 
ordinary Fokker-Planck equation (^) with the heuristic argument that the 
fluctuations of the gravitational potential during violent relaxation play the 
same role as collisions. With the additional (heuristic) argument that D ~ 
for large velocities, he could obtain a stationary solution of a Michie-King 
type §: 

- ( A(e-^-e-^ )e < emj 

f -\ e>e m . {bb) 

Since stars with e > e m are removed by the tidal field, this distribution 
function provides a depletion of high energy states and solves the infinite 
mass problem. 

Our present approach justifies the two phenomenological arguments of 
Lynden-Bell since equation (|54|) reduces to the Fokker-Planck equation if we 
assume / <C 770 (no degeneracy) and equation ([55]) gives a diffusion coefficient 
D ~ for large velocities. There is however a kind of "gap" in Lynden-Bell's 
approach since equation ( |56| ) does not reduce to the Fermi-Dirac statistics 
L6[) in the limit of low energies. Using the more general equation (B4) which 
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properly accounts for degeneracy effects, we now try to build up a "one piece" 
distribution function which makes the bridge between Lynden-Bell's statistics 
( p^| ) for e <C e m and the Michie-King model ( |56| ) for e ~ e m . 

During the stage of violent relaxation, the stars extract their energy from 
the rapid fluctuations of the gravitational field. By this process, some stars 
may acquire very high energies and escape from the system (being ultimately 
captured by the gravity of other systems). For these stars, D = ^ is a good 
approximation and equation (154) reads 



Of df -df d 
— + v— + F — = — 
dt dr <9v dv 



(57) 



We seek a stationary solution of the form / = /(e) 
^(^■) = (valid for sufficiently large |v|), we obtain 



Using the identity 



°-s 



K[% + Pttm-f) 



(58) 



or, equivalcntly, 



df 
de 



+ Pf(m - /) = -J, 



(59) 



where J is a constant of integration. If J = 0, we recover Lynden-Bell's 
distribution function (|l6|). If J ^ 0, equation (|59"| ) accounts physically for an 
escape of stars at a constant rate J. The system is therefore not truly static 
since it looses gradually stars but we may consider that the galaxy passes by 
a succession of quasi-stationary states which are solution of equation (^9|) . As 
stated previously, this equation is valid only for high energy stars. For lower 
energies, the system has settled down in a pure equilibrium state and J = 0, 
leading to Lynden-Bell's distribution. 

Our goal, now, is to solve equation (p9|). Put under the form 



df 
de 



+ Pvof~(3f +J = 0, 



(60) 



1 u' 



we recognize a Riccatti equation |33| . With the change of variables / — — ^ — , 
it can be converted into a linear partial differential equation of second order 



de 2 



,du 
~dl 



J(3u = 0. 



(61) 



The associated characteristic polynomial x 2 + t\q(5x — J (5 = has a strictly 
positive discriminant A = rfafi 2 + AJ(3 > 0. Therefore, the general solution of 
equation (|6l]) is 

— '-(Ae* e + Be~i e ), (62) 
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where 5 = v A and A, B are constants of integration. The solution of equation 
( p(j| ) is therefore 

1 or d — ~ e ■ ( 63 ) 

Setting A = A/B (B ^ otherwise / would be constant), it can be rewritten 



— _ 1 XjyoP - 5) + (r io/3 + 5)e- & 
2(3 A + e 



This distribution function vanishes at the escape energy e — e m defined by 

A(?7o/3 - 5) + (770/3 + 8) e - 5t ™ - 0. (65) 
With this new variable, we obtain the result 

7=Ay?0 (A- e e -^)(A+"e^) ' (66) 
where 5 is a solution of 

S = ^ x-l~-^ - (67) 

Now, for the cases of physical interest A ^> e~ /3l,oem , which means that de- 
generacy is negligible for energies close to the escape energy e m . In that case, 
S ~ rja/3 and we obtain the final result [ O : 

_ g-/3i?oc _ e -/3>)oe m 

f = ??0 : s • (68) 

The previous calculation is justified as long as D ^ 0, corresponding to 
relatively strong fluctuations. When the fluctuations die down at the end of 
the relaxation, the diffusivity D and therefore the diffusion current |J| go to 
zero. There is no more evaporation but the distribution function ( |68|) is main- 
tained as a stationary solution of the Vlasov equation. When e ~ e m , we re- 
cover the Michie-King model (|5rj| ) and when e<e m equation ([38]) reduces to 
the Fermi-Dirac distribution function (|f6|). Therefore, the distribution func- 
tion ( p8| ) connects continuously the two limits considered by Lynden-Bell [ [43| 
and can serve as a relevant model for (possibly degenerate) collisionless stel- 
lar systems. In particular, this distribution function could describe galactic 
halos limited in extension as a consequence of tidal interactions with other 



systems 34 . It could also be used to calculate realistic equilibrium states of 



collisionless stellar systems without the artifice of a material box. However, 



the main results of the box model 1 25 should not be dramatically altered 



The previous model describes self-gravitating systems subject to tidal 
forces. This form of confinement was first introduced in the case of globular 
clusters tidally trucated by a nearby galaxy [B[ . This model can also describe 
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a fraction of elliptical galaxies living in a rich environement. However, for the 
majority of elliptical galaxies, tidal forces are weak and the galaxy can be 
considered as isolated. Now, for isolated systems, other processes can account 
for "incomplete relaxation" and lead to different truncated models. Starting 
from equation (^), we can show that the diffusion coefficient is proportional 
to the fluctuations of the distribution function integrated over the velocity: 



Since these fluctuations rapidly decay as we depart from the relatively well- 
mixed central region of the galaxy, the diffusion current decreases also and 
this results in a confinement of the distribution function. It is expected there- 
fore that the relaxation will be effective only in a limited region of space where 
f 2 — / is sufficiently large. On the other hand, as the system develops finer 
and finer filaments during the mixing process, the diffusion coefficient is ex- 
pected to decrease in time. The diffusion current takes therefore increasingly 
small values and the relaxation is slowed down. This decay may be quite 
rapid and a quantitative treatment of this effect would require a better un- 
derstanding of the correlation function /(r, v, i)/(r', v', t) whose expression 
was simply postulated in section ||. The development of these ideas will lead 
to other truncated models probably closer to those introduced phenomeno- 
logically by, e.g., Stiavelli & Bertin |Q, Tremaine |^9| and Hjorth & Madsen 
| |3l) . In particular, the above discussion is quite consistent with the scenario 
of incomplete violent relaxation developed by Hjorth & Madsen. These au- 
thors introduce a two-step process: (i) in a first step, they assume that violent 
relaxation proceeds to completion in a finite spatial region, of radius r max , 
which represents roughly the core of the galaxy (where the fluctuations are 
important). At this stage, the escape energy represents no special threshold 
so that negative as well as positive energy states are populated in that re- 
gion, (ii) After the relaxation process is over, positive energy particles leave 
the system and particles with ^(r maa; ) < e < move in orbits beyond r max , 
thereby changing the distribution function to something significantly 'thin- 
ner' than a Boltzmann distribution. The crucial point to realize is that the 
differential energy distribution N(e), where N(e)de is the number of stars 
with energy between e and e + de, will be discontinuous at the escape energy 
e = since there is a finite number of particles within r max after the relax- 
ation process. It can be shown that this discontinuity implies necessarily that 
/(e) ~ (— e) 5 / 2 for e — * 0~ p5| . Therefore, the truncated model consistent 
with this scenario is [Blfl: 




(69) 




(70) 



This truncated model corresponds formally to a composite configuration 
made of an isothermal core and a polytropic envelope with index n = 4. 
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Of course, if the core is degenerate, the Boltzmann distribution for e < e' 
must be replaced by the Fermi-Dirac one. This truncated model gives very 
good fit with elliptical galaxies and can reproduce the i? 1 / 4 law |3^]. The box 
model of section |] can be considered as a simple approximation of this more 
realistic model, the "box" playing the role of the envelope. In fact, the sce- 
nario developed by Hjorth & Madsen is almost equivalent to complete violent 
relaxation in a finite container with the container removed after the relax- 
ation. This scenario solves the infinite mass problem and rehabilitates the 
use of statistical mechanics to understand the structure of elliptical galaxies. 



8 Conclusion 

We have described in this paper the process of violent relaxation in stellar sys- 
tems from the viewpoint of statistical mechanics, as originally introduced by 
Lynden-Bell |Q . Lynden-Bell proposed that the structure of galaxies could 
result from a law of chaos: there is a total lack of information at small scales 
since the stars have complicated orbits, but the exciting phenomenon is that 
microscopic disorder leads to macroscopic order. The same ideas of statisti- 
cal mechanics have been introduced in two-dimensional turbulence described 
by the Euler equation p3p^p7| , p3t to explain the formation and mainte- 
nance of large scale coherent vortices like the Great Red Spot of Jupiter or 
the cyclones and anticyclones that populate the earth atmosphere |2^,^,^|. 
The formation of self-organized vortices in two-dimensional turbulence can 
also have applications in the context of planet formation where large-scale 
vortices present in the solar nebula could efficiently trap dust particles to 
form the planetesimals and the planets (see a complete list of references in 
Chavanis In fact, the statistical mechanics of the 2D Euler equation is 

equivalent to the theory of Lynden-Bell applying to the Vlasov equation. In 
a sense, the Vlasov equation is just the Euler equation for a "fluid" evolving 
in a six-dimensional phase space. In this analogy Jll]|l4|,[l7]] , the vorticity and 
the stream function in 2D turbulence play the same role as the distribution 
function and the gravitational field in galaxies. This analogy concerns not 
only the equilibrium states (the formation of large-scale structures) but also 
the relaxation towards equilibrium ]27j and the statistics of fluctuations |24| . 
A kinetic theory of two-dimensional vortices can be developed in analogy 
with stellar dynamics [ fl9|| . In this kinetic theory, a point vortex experiences 
a diffusion process and a "systematic drift". This "systematic drift" is 
the counterpart of the "dynamical friction" |H][ ] experienced by a star as a 
result of close encounters. Relaxation equations analogous to equations (Blj) 
and (|5^) have been derived in the context of vortex dynamics [ p^T^Jl^ , ^9{ . 
In addition, the problem of "incomplete relaxation" is also encountered in 2D 
turbulence to explain the confinement of a vortex (e.g., a dipole or a tripole) 
that forms after a rapid merging |27],^3],[32| . It has been demonstrated ex- 
plicitly in two-dimensional turbulence (where the numerical simulations are 
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easier to implement) that the relaxation equations derived from the MEPP 
and including a space dependant diffusion coefficient related to the fluctu- 
ations of the vorticity (analogous to Eq. (69)) can account for this kinetic 
confinement [^2|. We believe that the relaxation equations proposed in the 
stellar context should work as well. The statistical mechanics of continuous 
vorticity fields also predicts a Fermi-Dirac distribution at equilibrium with the 
same interpretation of the degeneracy as in Lyndcn-BclPs theory. Although 
this degeneracy is hard to evidence for galaxies (and remains controversal) , it 
has been vindicated by various numerical simulations and laboratory exper- 
iments of two-dimensional fluids. This suggests that the degenerate version 
of the theory should also be used in the stellar context, ff all these effects 
are taken into account correctly, it is plausible that the statistical mechanics 
of 2D vortices and self-gravitating systems has a chance to account for the 
fascinating process of self-organization in nature. 
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